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SUMMARY 

This paper describes the implementation of an improved quadrilateral 
membrane element in NASTRAN, Descriptions of the geometrical and kinematic 
properties of the element are included along with development of the matrices 
and vectors which characterize the element. The necessity for considering 
small deviations from planeness of the element is discussed and the approach 
taken to account for these deviations for the new element is described. The 
improved accuracy available from the element over the existing quadrilateral 
membrane is indicated by a sample calculation for which an analytical solution 
is available from beam theory. For the same finite element idealization, the 
errors in maximum displacement and stress were significantly reduced by the 
use of the new element. 


INTRODUCTION 

One of the more frequently used elements in the NASTRAN library is the 
quadrilateral membrane element (CQPMEM) . This element is used to represent 
portions of structures for which membrane action constitutes the predominant 
contribution to the strain energy. An example of one such application is to 
the skin of aircraft wings. In addition, the quadrilateral membrane element 
is also used in conjunction with the quadrilateral bending element (CQDPLT) to 
form the membrane-bending elements CQUAD1 and CQUAD2 needed to represent more 
general deformation behavior. 

The quadrilateral membrane element currently available in NASTRAN 
is composed of overlapping constant strain triangles as described in reference 
1, It has been reported that the CQDMEM element does not accurately represent 
problems involving high stress gradients, suggesting that the current element- 
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needs improvement. The reason for this difficulty is generally attributed to 
the constant strain field provided by this element (ref. 2). Several improved 
quadrilateral membrane elements have been developed which have linear strain 
fields (e.g. refs. 3 and 4) and these elements should provide improved results 
in membrane element applications. The implementation of an improved quadri- 
lateral membrane element in NASTRAN was undertaken as an in-house project to 
improve the element library by attempting to overcome the noted shortcomings. 

The element chosen for implementation is the linear strain isoparametric 
quadrilateral membrane element described in references 3 and k. The reasons 
for this choice are: 

(1) The element is conforming, i.e. the displacements of adjoining 
elements are matched along their entire interface. As a result, 
the strain energy is an upper bound to the corresponding exact 
energy . 

(2) The stresses and strains vary within the element thus providing 
improved accuracy over the existing element. 

(3) The element is well-documented. 

The purposes of the present report are to present a description of the 
element being implemented and to demonstrate the increased accuracy available 
through the new element by comparison with the existing quadrilateral element 
and an analytical solution. 


SYMBOLS 


[A] 

[B] 

Ec] 

[E] 


matrix relating strains and displacements 
transformation matrix relating displacements in mean plane 
to those at actual grid points 

matrix which relates rotations of the mean plane to displacements 
of mean plane 

matrix relating displacements in element coordinate system 
to those in basic coordinate system 
membrane strains 

vector of forces at actual grid points in element coordinate 
system 
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<f e } 


[G e ] 


H 

h 

J 


[K] 

[K^] 


[K ee^ 

[K^ e ] 


{P} 

{ V 

[S] 

[s e l 

{S t } 


[T] 


u,v,w 

u,v,w 




V 

X.Y.Z 

x,y,z 


vector of forces in mean plane 
matrix relating stresses and strains 
distance from actual grid point to the mean plane 
thickness of membrane 

Jacobian of transformation from x-y coordinates to 
coordinates 

stiffness matrix referred to global coordinate system 
differential stiffness matrix referred to global coordinate 
system 

stiffness matrix in element coordinate system 
differential stiffness matrix in element coordinate system 
thermal load vector referred to global coordinate system 
thermal load vector in element coordinate system 
matrix relating element stresses to global displacements 
matrix relating element stresses to element displacements 
vector relating element temperature to element stress 
basic to global coordinate transformation matrix 
reference or stress-free temperature of the element 
temperature of the element above the reference or stress-free 
temperature 

displacements in x-, y-, and z-directions, respectively 
displacements in X-, Y-, and Z-directions, respectively, see 
table I 

vector of displacements at actual, grid point in element 
coordinate system 

vector of displacements in mean plane 

cartesian coordinates used in table I 

element cartesian coordinate directions, see figure 1 

element parametric coordinates, see figure 1 

rotations of the mean plane element about x,y and z axes, 

respectively 


Subscripts: 

1* 2, 3, 4 refer to grid points 1, 2, 3 and 4 respectively, of the element 
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A subscript preceded by a comma indicates partial differentiation with 
respect to the subscript. 


DESCRIPTION OF THE ELEMENT 


In this section of the paper, descriptions of the geometry and kinematic 
behavior of the isoparametric element (figure l) will be given. The element, 
since it is defined by four points, need not be planar j however the development 
of the necessary matrices is carried out for a flat element. The treatment of 
the case for which the four points are not coplanar will be discussed in a 
later section of the paper. The element parametric coodinates £ and r) 
shown in figure 1 vary linearly between zero and one with the extreme values 
occurring on the sides of the quadrilateral. Further, lines of constant £ 
and n are straight as indicated in the figure. A set of element cartesian 
coordinates x, y, z is defined as follows: the x-axis is along the line 

connecting the first two grid points; the y-axis is perpendicular to the 
x-axis and lies in the plane of the element; and the z-axis is normal to the 
plane and forms a right handed system with the x- and y-axes. Displacement 
components in x, y, and z directions are denoted by u, v, and w, 
respectively. As given in reference 3, the displacement field is assumed to 
have the following form: 


u(£,n) ■ (1 - 0(1 - n) u x + £(i - n) u 2 + £nu 3 + (1 - 5) rp^ 


v(£,n) ■ (1 - Od - n) + ?(i - n) v 2 + £nv 3 + (l - O nv^ 


(i) 


where the subscript on a displacement component denotes the grid point value 
of the component. 

It may be observed that on lines of constant £, u and v vary linearly 
with r\ and on lines of constant ri, u and v vary linearly with £. In 
particular u and v vary linearly on each edge between grid points and as 
a result the displacements of adjacent elements are matched all along their 
common edges. The element is therefore of the "conforming" type which 
guarantees that the element will converge as an upper bound on the strain energy. 
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The required membrane strains are related to the displacements u and v 
by the familiar relations 


e 


x 


u. 



e 

xy 



( 2 ) 


where the comma indicates partial differentiation. Since the displacements 
are expressed in terms of £ and n» the operations in equations (2) cannot 
be carried out without knowing the relationships between the (x,y) coordinates 
and the (£,ri) coordinates. These relations are given in reference 3 as 


x = (i - 5)(i - rO x x + s(i - n) x 2 + Cnx 3 + (i - Z) 
y - (i - U(i - n) y x + 5d - n) y 2 + Sny 3 + (1 - O ny^ 


(3) 


By use of familiar relations involving partial derivatives, the operations 
indicated in equations (2) may be performed. Thus, for example, 


u 


’x 


u, £ ?, x + u *n n, x 


(U) 


where 


z 


»x 



n.. 



(5) 


and 


J = 




( 6 ) 


P’E y 'n| 

It is noted in passing that equations (l) which relate displacements within 
the element to its grid point values or "parameters" are identical in form 
to equations (3) for coordinates x and y. Thus the term "isoparametric" 
is applied to characterize the element. 

For the special case of a rectangle it can be shown that the x and £ 
directions are identical as are the y and n directions. In this case e x 
is linear with respect to y and constant with respect to x, whereas e 

y 

is linear with respect to x and constant with respect to y. The shear 
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strain varies linearly with respect to hoth x and y. In contrast to the 
strains, all three stress components vary linearly in hoth the x and y 
directions. This is a direct consequence of the constitutive equations. For 
nonrect angular shapes the behavior of the stress and strain components is 
more complicated and is not easily characterized. 

IMPLEMENTATION OF THE ELEMENT 

The addition of an element to NASTRAN requires the derivation of a set 
of characteristic matrices and vectors as described in reference (l) and those 
necessary for the new element axe: 

stiffness matrix, [K ] 

lumped mass matrix, [M ge ] 

thermal load vector, {P } 

e 

stress recovery matrices, [S ] and {S^} 

e t 

differential stiffness matrix, [Kr ] 

The development of these matrices is presented in this section along with a 
description of the procedure used when the four grid points defining the 
element are not coplanar. Finally, this section contains an outline of the 
matrix transformations required so that the new element will be compatible 
with others in NASTRAN. 

Stiffness Matrix 

Using equations (l) through (5) results in the following relation between 
strains and grid point displacements 
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where the elements of the 3x8 matrix [A] are functions of £ and n* The 
stress-strain relation is given by 



where the 3x3 matrix [G ] represents most generally a completely anisotropic 
material. The terms a x , and a_^ are thermal expansion coefficients 
and t is the average temperature of the element above the stress-free 
temperature, given by 

* ■ r <*i + % * s * V - *0 (9) 

The strain energy, V, (apart from thermal effects) is 

v - | /q Jq <»/ [A] T [G e ][A]{u e } J d? dn (10) 

and from this expression, the stiffness matrix can be identified as 

[K ee I » h /J [A] T [G e ][A] J d5 dr) (11) 
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The required integration is performed numerically by the use of Gaussian 
quadrature using a 4x4 grid. For a discussion of Gaussian quadrature as used 
for isoparametric elements, see references 3 and 4. 

Lumped Mass Matrix 

The mass matrix developed for the isoparametric element is a lumped mass 
matrix since coupled mass matrices for membrane elements generally result in 
overly stiff representations in dynamic problems (see ref. 1 p, 5.5-5). One 
method of lumping is to assign one quarter of the mass of the element to each 
of the four grid points. However , this method usually does not preserve 
the location of the center of mass of the element. Accordingly, the method 
used to generate the mass matrix for the new element is that presently used 
in NASTRAN for the existing quadrilateral membrane element. This later method 
is based on an averaging procedure which always preserves the center of mass. 


Thermal Load Vector 


For the purpose of developing the thermal load vector, the contribution 
to the potential energy, U, of the element .temperature above some stress-free 
value is written as 


U 


= h M 


r ^ 

T 

r 

a 

a I 

X 


X 

a 

L < 

a l 

y 


y 

a 


a 

xy 


xy 

V J 


J 


t j dUn 


( 12 ) 


Using equations (7) and (8) in equation (12) and discarding an irrelevant 
constant term not involved in the solution, results in 


» - » *1 'o ‘V T [A]V e ] 


r a 

X 


<v 


t J d£ dr) 


(13) 


and the thermal load vector is then recognized as 


a 

xy 
V. J 
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r 


{P e> ’ h 4 4 U)T 


a 


J d5 dn Ca e 3 / a y 


(1U) 


Stress Recovery Matrices 

Expressions for the element stress components o 


and a at any 


x' y ' xy 

point written in terms of displacements measured in the element cartesian 
coordinate system are obtained by combining equations (7) and (8) to give 


where 



tW - to > 1 


( 15 ) 


r 

a x 

[S e ] * [G- e ] [A] and {S^.} * [G g ] < a y \ (l6) 

a 

L*yJ 

Although the stress components may vary within the element, they are computed 
only at the intersection of the element diagonals for the purpose of stress 
output. Once the stress components are obtained, they are used to compute the 
principal stresses and directions by appropriate formulas given in reference 5, 


Differential Stiffness Matrix 

The differential stiffness matrix, which is used in NASTRAI primarily 
for linear buckling analyses, is developed by a consideration of the work 
done by stress components ct x , a y , and during small rotations co x , 

w y , and w z about the three element cartesian axes. The expression for the 
work done is taken from reference 1 (section 7.1, equation (16)) given as 
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w = - 


h <*1 |*1 
2 J 0 J 0 


r 

0) 

x 

T 


f ^ 
0) 

X 

03 

y 

> 

[K d ]) 
m N 

0) 

y 

1 

L j 



03 

Z 

V. J 


(17) 


where 


[K d ] 

cow 


-0 


xy 


xy 


0 0 

and the rotation components are given by 


0 

0 


a +0 

x y_ 


(18) 


C0„ = W, CO = -w, 
x y y x 


“z * I <v 'x - “’y 1 


(19) 


In order to evaluate co and co , the behavior of w in the element is 

x y* 

required, and for this purpose w is assumed to have the same parametric 
variation as u and v, thus 


v » (i - 5)(i - n) + c(i - n) w £ + 5nw 3 + (l - £) 


Combining equations (l), (19 ) , and (20) results in 


( 20 ) 


r 

co 


(to 


} = [C]{u e > 


( 21 ) 


W z 


where the elements of the 3x8 matrix [C] are functions of £ and fj. 
Substituting equation (21) into equation (17) gives 
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( 22 ) 


« - - | 'o 4 {u e )T t C l T lO tC]tu e> J « 

and the differential stiffness matrix is recognized as 

[K Je ] = h f 0 f l [C]T [K wJ [C] J dC dn (23) 

The integration is performed by use of Gaussian quadrature using a 4x4 grid. 

The values of cr x , , and a at each of the sixteen quadrature points are 
obtained from equation (15 ) for appropriate values of £ and n. 

Use of the Mean Plane 

As mentioned previously, the matrices and vectors associated with the new 
quadrilateral membrane element were derived for a planar element. If the four 
grid points are not coplanar they are projected onto a so-called mean plane 
which is defined in the following manner (see figure 2), Two skewed lines 
whose end points are grid points 1 and 3 and grid points 2 and 4, respectively, 
are defined to be the diagonals of the element. The mean plane is defined 
so that it passes through the midpoint of the perpendicular connector of the 
two diagonals and parallel to them (see ref. 5 p. 4.87-105.). As shown in 
figure 2, if the length of the perpendicular connector is 2H, then the user 
defined grid points are alternatively H units above and below the mean plane 
as one progressively moves around the element. Once the actual grid points 
are projected onto a mean plane, a planar element is defined and the previously 
derived matrices are applicable. 

The authors would like to emphasize that the use of the mean plane concept 
does not finally resolve the question of how to deal with a quadrilateral 
element whose grid points are not coplanar and that further research on 
this subject is warranted. One alternative to the mean plane is to 
assume the element to be planar with the plane defined by three of the grid 
points. This alternative was found to be undesirable for the new element 
because numerical results were sensitive to the grid point numbering sequence 
in a single element. This numerical experiment to determine this sensitivity 
is described in the Appendix. 
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Transformation of Matrices and Vectors 


to Global Coordinates 

Up to this point the characteristic matrices and vectors for the element 
have been derived in terms of displacements in the element cartesian 
coordinate system. In order that the new element be compatible with other 
elements in NASTRAN , it is necessary to transform the matrices and vectors 
so that they are expressed in the global coordinate system. This procedure 
involves three transformations: 

(1) transformation from displacements in the mean plane in the element 
cartesian system to the displacements at the user-defined grid points 
in the element cartesian system 

(2) transformation from displacements in the element coordinate system to 
displacements in the NASTRAN basic coordinate system 

(3) transformation from the NASTRAN basic coordinate system to the NASTRAN 
global coordinate system. 

The first transformation is based on replacing the set of forces in the 

mean plane {f } by a statically equivalent set at the user defined grid points 
6 

{f }. The first set consists of forces in. the x and y directions only, 
whereas the second set consists of forces in the x, y, and z directions. 

This statement of static equivalence can be written as 


where 


{f } = [B]{f } 
a e 


r ^ 

f xl 


V| 

f yi 


V 

f zl 


f x2 

W H> 

TO 

and {f } = 

f y2 

• 


f x3 

• 


f y3 

• 


^xh 

f zU 

v. _ 

/ 

V 

V 


( 2*0 


(25) 
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The matrix [B] depends on the geometry of the element and H which is a 
measure of the nonplaneness of the element. This same matrix is used to 
transform displacements, thus 


{u > = [B] T {u } (26) 

“ cl 

The second and third transformations are required of ail elements in NASTRAN 
and are discussed fully in reference 1. In terms of the notation of reference 
1 the transformation matrix from element coordinates to NASTRAN basic 
coordinates is denoted by a 12x12 matrix [E] , and the transformation matrix 
from NASTRAN basic coordinates to NASTRAN global coordinates is denoted by a 
12x12 matrix [T], The relation between element and global displacements may 
finally be written as 


{u } = [B] T [E] T [T] {u } (27) 

e g 

All that is now required to express the matrices and vectors in the global 
coordinate system is to substitute equation (27) into equations (10), (13), 
(15), and (22). The results are summarized below 

Stiffness matrix : 


[K] = [T] T [E][B][K ee ][B] T [E] T [T] (28) 


Thermal load vector: 


{P} = [T] T [E][B]{P e ) (29) 

Stress recovery matrix : 

[s] = [S ][B] T [E] T [T] ( 30 ) 

e 

Differential stiffness matrix; 


(31) 


[K d ] = [T] T [E][B][K d ][B] T [E] T [T] 
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PRELIMINARY RESULTS 


At the time of writing, the element is in the early stage of being veri- 
fied. A stand-alone program has been written for the purpose of testing the 
element in a NASTRAN-like environment . The first available test case for this 
element is a cantilever beam modeled by 1.6 equal elements as shown in figure 3. 
Displacements and stresses were computed using the existing CQDMEM element and 
the new isoparametric element based on the same finite element idealization. 

The finite element results are then compared to an analytical solution based 
on elementary beam theory. In figure 3 the analytical solution for displace- 
ments is indicated by the solid line and the finite element results by dashed 
lines . The new element is seen to be a significant improvement over the CQDMEM 
element. When comparison is made with the analytical solution, the largest 
error exhibited by CQDMEM is 50% whereas the corresponding error for the new 
element is 9%. The comparison of stresses is presented in figure 1+ where 
the results have been normalized to the maximum stress predicted by elementary 
beam theory. Although the finite element stresses were computed at a single 
point in each element, these stresses are presented as continuous curves. 

Again the exact solution is shown by the solid line and the finite element 
results by dashed lines. The results indicate a significant improvement in 
accuracy by using the new element. The maximum error in stress using the 
CQDMEM element is about 1+9% whereas that for the improved element is lb%, 

CONCLUDING REMARKS 

An isoparametric quadrilateral membrane element which is being 
implemented in NASTRAN by the NASTRAN Systems Management Office is described. 
Included are descriptions of the basic geometry and deformation characteristics 
of the element, as well as an outline of the derivation of each of the required 
matrices and vectors for the element. Derivations of the stiffness matrix, 
thermal load vector, stress recovery matrices, and differential stiffness 
matrices are given. A discussion of the need for some method of accounting 
for deviations from planeness is included along with a description of the 
method used for this purpose in the new element. Finally, a sample 
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calculation is carried out in which the displacements and stresses in a 
cantilever beam are computed by the new element and by the existing quadri- 
lateral membrane (CQDMEM) and compared with an analytical solution from beam 
theory. The significant increase in accuracy obtainable with the new element 
is demonstrated by the calculation, 
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APPENDIX 


SENSITIVITY OF MEMBRANE RESULTS TO 
ELEMENT GRID POINT NUMBERING SEQUENCE 

A numerical experiment suggested "by Dr. Raphael Haftka, currently a 
National Research Council postdoctoral Fellow at the Langley Research Center, 
is performed on the structure sketched at the top of table I. The top 
member is non-planar and the vertical member is planar. Both members were 
represented by the new quadrilateral membrane element assumed to be planar 
with the plane of each element being defined by three of its grid points. 
Displacements under two loading conditions were computed by means of a 
stand-alone program. In the calculations, the grid point sequence for the 
top member was varied as shown in table I, and the effects of these varia- 
tions are shown in table II. Examination of table II (a) shows that for a 
load applied in the X-direction, the displacements both in the X-Y plane and 
normal to it are sensitive to the grid point numbering. For example, there 
is a 12% difference among the displacements in the Y-direction as well as 
among the displacements in the Z-di recti on at grid point 3. For the case of 
loads in the Z-direction there is a difference among the displacements in 
the Y-direction of 46%. The consequences of these percentage errors are 
tempered by the observation that the large errors noted always occurred in 
displacement components which were not in the direction of the applied load. 
By contrast the maximum difference was 1.5% for a displacement in the 
X-direction at grid point 2 in the first loading condition. Although the 
difference for displacements in the direction of the load is small, it is 
felt that the existence of a situation where different answers can be 
obtained with the same finite element model by merely changing the grid 
point numbering sequence for individual elements should be avoided. 
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TABLE I.- STRUCTURE USED TO DETERMINE EFFECT OF NONPLANENESS 



(a) Gridpoint Locations 
(all coordinates given in centimeters) 


Gridpoint 

X 

Y 

z 

1 

1280.45 

480.44 

9.38 

2 

1367.30 

480.44 

7.36 

3 

1391.87 

536.68 

5.52 

4 

1305. OX 

536.68 

5.78 

5 

1367.30 

480.44 

-87.00 

6 

1391,87 

536.68 

-87.00 


(b) Loading 


Loading condition 1 

F X2= F X3- 1N 

Loading condition 2 

F 72 = F Z3 = 1N 


(e) Grid Numbering Sequence for Top Element 


Case 

1st Point 

2nd Point 

3rd Point 

4th Point 

I 

1 

2 

3 

4 

II 

2 

3 

4 

1 

III 

3 

4 

1 

2 

IV 

4 

1 

2 

3 
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Table XI.- Effect of Grid Point Numbering on Displacement of 

a Warped Quadrilateral. (All Displacements are in cm.) 


(a) Displacement due to loads in X-direction 

Grid point 
Sequence (table I) 

\ x 10 6 

v 2 x 10^ 

v 2 x 10 6 

x 10 

X 10^ 

— 6 
x 10 

I 

5.46 

1.59 

1.1+6 

6.85 

.1+98 

-1.07 

II 

5.1+9 

1.57 

1.1+8 

6.81+ 

.526 

-1.08- 

III 

5.1+1+ 

1.56 

1.35 

6.86 

• 511 

-1.20 

IV 

5.1+1 

1.57 

1.33 

6.88 

.1+70 

-1.19 


(b) Displacements due to loads in Z-direction 

Grid point 
Sequence 

u 2 x 10 6 

Vg x 10 ^ 

Wg x 10 ^ 

u^ x 10 ^ 

v 3 Xl ° 

W 3 x 10 ^ 

I 

.1+85 

• 535 

6 . 1+1 

-.095 

-.298 

6.05 

II 

.506 

.365 

6.33 

-.099 

-.1+1+7 

6.12 

III 

.377 

. 1+21 

6.33 

-.223 

-.392 

6.12 

IV 

.351+ 

.591 

6 . 1+1 

-.217 

-. 21+6 

6.05 
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Mean p. 
points 
grid pi 







Young’s Modulus * 100 GN/m^ 
Poisson’s ratio = 1/3 


— exact solution 

— isoparametric element 

— existing GQDMEM element 








